# Clean workspace
rm(list = ls())  # Remove all objects
graphics.off()    # Close all graphics devices
cat("\014")      # Clear console

# Install required packages if not already installed
if (!require(tidyverse)) install.packages("tidyverse")
if (!require(readxl)) install.packages("readxl")
if (!require(lme4)) install.packages("lme4")
if (!require(moments)) install.packages("moments")
if (!require(writexl)) install.packages("writexl")  # Add writexl package for Excel output

# Load required libraries
library(tidyverse)  # for data manipulation
library(readxl)     # for reading Excel files
library(lme4)       # for mixed effects models
library(moments)    # for skewness and kurtosis
library(writexl)    # for writing Excel files

# Function to load data
load_erp_data <- function(file_path, sheet_number = 3) {  # Default to load the 3rd sheet
  # Determine file extension
  file_ext <- tools::file_ext(file_path)
  
  # Load data based on file type
  data <- if (file_ext == "csv") {
    read.csv(file_path)
  } else if (file_ext %in% c("xlsx", "xls")) {
    read_excel(file_path, sheet = sheet_number)  # Specify sheet
  } else {
    stop("Unsupported file format. Please provide a CSV or Excel file.")
  }
  
  # Convert to tibble for better handling
  data <- as_tibble(data)
  
  # Ensure proper column types
  data <- data %>%
    mutate(
      Subject = as.factor(Subject),
      Semantic_Relatedness = as.factor(`Semantic Relatedness`),
      Priming_Word_Type = as.factor(`Priming Word Type`),
      Channel = as.factor(Channel),
      Group = as.factor(Group)
    )
  
  return(data)
}

# Set input file path
input_file <- "source_data/Ch_Ja_Lex_ERP_2Group_for R.xlsx"

# Create result folder name from input file name
result_folder <- tools::file_path_sans_ext(basename(input_file))
result_folder <- gsub("_for R", "", result_folder)  # Remove "_for R" suffix
result_path <- file.path("results", result_folder)

# Create result folder
if (!dir.exists(result_path)) {
  dir.create(result_path, recursive = TRUE)
  print(paste("Created result folder:", result_path))
}

# Load your data (3rd sheet)
print("Loading data...")
erp_data <- load_erp_data(input_file, sheet_number = 3)
print("Data loading completed")

# View data structure
print("\n============= Data Structure =============")
print("Data dimensions:")
print(dim(erp_data))
print("\nVariable names:")
print(names(erp_data))
print("\nData structure:")
str(erp_data)

# View ERP time windows
print("\n============= ERP Time Windows =============")
# Get all column names starting with "ERP_"
erp_columns <- names(erp_data)[grep("^ERP_", names(erp_data))]
print("ERP time windows list:")
print(erp_columns)

# Create a data frame to store all results
results_df <- data.frame(
  Time_Window = character(),
  Mean = numeric(),
  SD = numeric(),
  Skewness = numeric(),
  Kurtosis = numeric(),
  Skewness_Acceptable = logical(),
  Kurtosis_Acceptable = logical(),
  Mean_Acceptable = logical(),
  Skewness_Criteria = character(),
  Kurtosis_Criteria = character(),
  Mean_Criteria = character(),
  Kurtosis_Interpretation = character(),
  Kurtosis_Detail = character(),
  Extreme_Values_Count = integer(),
  Extreme_Values_Detail = character(),
  Cleaned_Data_Kurtosis = numeric(),
  Cleaned_Data_N = integer(),
  Removed_Data_Percentage = numeric(),
  stringsAsFactors = FALSE
)

# View basic statistics for each time window
print("\nBasic statistics for each time window:")
for (col in erp_columns) {
  print(paste("\n", col, ":"))
  print(summary(erp_data[[col]]))
}

# View factor levels
print("\n============= Factor Levels =============")
print("Number of subjects:")
print(length(unique(erp_data$Subject)))
print("\nSemantic Relatedness levels:")
print(levels(erp_data$Semantic_Relatedness))
print("\nPriming Word Type levels:")
print(levels(erp_data$Priming_Word_Type))
print("\nChannels:")
print(levels(erp_data$Channel))
print("\nGroups:")
print(levels(erp_data$Group))

# Function to get detailed interpretation for normality
get_normality_interpretation <- function(skewness, kurtosis, mean_val) {
  skew_text <- if(abs(skewness) < 1) {
    "偏度正常，分布对称性良好"
  } else if(skewness > 1) {
    "正偏度，分布右侧有长尾"
  } else {
    "负偏度，分布左侧有长尾"
  }
  
  kurt_text <- if(abs(kurtosis) < 2) {
    "峰度正常，分布形态接近正态"
  } else if(kurtosis > 2) {
    "峰度过高，分布较尖峭（重尾），极端值较多"
  } else {
    "峰度过低，分布较平坦（轻尾），极端值较少"
  }
  
  mean_text <- if(abs(mean_val) < 0.1) {
    "残差均值接近0，无系统性偏差"
  } else {
    "残差均值偏离0，可能存在系统性偏差"
  }
  
  return(paste(skew_text, kurt_text, mean_text, sep = "; "))
}

# Model assumption tests for each ERP time window
print("\n============= Linear Mixed Effects Model Assumption Tests =============")

# Create a function to check model assumptions
check_model_assumptions <- function(model, time_window) {
  residuals <- resid(model)
  
  # 获取原始数据
  original_data <- data.frame(
    Time_Window = time_window,
    Residual = residuals,
    Subject = erp_data$Subject,
    Channel = erp_data$Channel,
    Group = erp_data$Group,
    Semantic_Relatedness = erp_data$Semantic_Relatedness,
    Priming_Word_Type = erp_data$Priming_Word_Type,
    ERP_Value = erp_data[[time_window]]
  )
  
  print(paste("\nModel assumption test results for:", time_window))
  print("\n1. Residual descriptive statistics")
  print("Sample size:")
  print(length(residuals))
  print("\nMean:")
  print(mean(residuals))
  print("Standard deviation:")
  print(sd(residuals))
  print("Skewness:")
  skew_val <- skewness(residuals)
  print(skew_val)
  print("Kurtosis:")
  kurt_val <- kurtosis(residuals)
  print(kurt_val)
  
  # 初始化变量
  current_data <- original_data
  all_extreme_data <- data.frame()
  iteration <- 1
  max_iterations <- 5  # 增加到5次迭代
  # 使用更严格的IQR倍数序列
  iqr_multipliers <- c(2.5, 2.0, 1.5, 1.2, 1.0)  # 逐渐降低IQR倍数，最后使用1.0作为最严格标准
  kurtosis_threshold <- 2  # 可接受的峰度阈值
  min_sample_size <- nrow(original_data) * 0.7  # 设置最小样本量为原始数据的70%
  
  # 迭代清理极端值
  while (abs(kurtosis(current_data$Residual)) >= kurtosis_threshold && 
         iteration <= max_iterations && 
         nrow(current_data) > min_sample_size) {  # 添加样本量检查
    # 使用当前迭代的IQR倍数
    current_multiplier <- iqr_multipliers[iteration]
    
    # 计算当前数据的四分位数和IQR
    q1 <- quantile(current_data$Residual, 0.25)
    q3 <- quantile(current_data$Residual, 0.75)
    iqr <- q3 - q1
    lower_bound <- q1 - current_multiplier * iqr
    upper_bound <- q3 + current_multiplier * iqr
    
    # 标记极端值
    is_extreme <- current_data$Residual < lower_bound | current_data$Residual > upper_bound
    
    if (sum(is_extreme) == 0) break  # 如果没有找到新的极端值，退出循环
    
    # 检查移除这些极端值后的样本量是否足够
    if (nrow(current_data) - sum(is_extreme) < min_sample_size) {
      print(sprintf("\n警告：移除极端值后样本量将低于最小要求（%.0f个），停止清理", min_sample_size))
      break
    }
    
    # 获取当前迭代的极端值
    current_extreme <- current_data[is_extreme, ]
    current_extreme$Iteration <- iteration
    current_extreme$IQR_Multiplier <- current_multiplier
    current_extreme$Z_Score <- (current_extreme$Residual - mean(current_data$Residual)) / sd(current_data$Residual)
    
    # 添加到所有极端值数据框中
    all_extreme_data <- rbind(all_extreme_data, current_extreme)
    
    # 更新当前数据集（移除极端值）
    current_data <- current_data[!is_extreme, ]
    
    # 计算并打印当前迭代的结果
    current_kurtosis <- kurtosis(current_data$Residual)
    print(sprintf("\n迭代 %d (IQR倍数: %.1f):", iteration, current_multiplier))
    print(sprintf("发现 %d 个极端值", nrow(current_extreme)))
    print(sprintf("清理后峰度: %.2f", current_kurtosis))
    print(sprintf("剩余样本量: %d (%.1f%%)", nrow(current_data), 
                 100 * nrow(current_data) / nrow(original_data)))
    
    iteration <- iteration + 1
  }
  
  # 准备极端值汇总
  if (nrow(all_extreme_data) > 0) {
    extreme_summary <- data.frame(
      Iteration = 1:max_iterations,
      IQR_Multiplier = iqr_multipliers,
      Extreme_Values_Found = sapply(1:max_iterations, function(i) sum(all_extreme_data$Iteration == i)),
      Kurtosis_After_Removal = sapply(1:max_iterations, function(i) {
        temp_data <- original_data[!original_data$Residual %in% all_extreme_data$Residual[all_extreme_data$Iteration <= i], ]
        if (nrow(temp_data) > 0) kurtosis(temp_data$Residual) else NA
      }),
      Remaining_Sample_Size = sapply(1:max_iterations, function(i) {
        temp_data <- original_data[!original_data$Residual %in% all_extreme_data$Residual[all_extreme_data$Iteration <= i], ]
        nrow(temp_data)
      })
    )
    
    # 添加每次迭代的极端值分布信息
    extreme_distribution <- data.frame(
      Iteration = 1:max_iterations,
      Positive_Extremes = sapply(1:max_iterations, function(i) 
        sum(all_extreme_data$Residual[all_extreme_data$Iteration == i] > 0, na.rm = TRUE)),
      Negative_Extremes = sapply(1:max_iterations, function(i) 
        sum(all_extreme_data$Residual[all_extreme_data$Iteration == i] < 0, na.rm = TRUE)),
      Max_Z_Score = sapply(1:max_iterations, function(i) 
        max(abs(all_extreme_data$Z_Score[all_extreme_data$Iteration == i]), na.rm = TRUE))
    )
    
    # 保存详细的极端值信息
    extreme_file <- file.path(result_path, paste0(time_window, "_kurtosis_extreme_values.xlsx"))
    write_xlsx(list(
      Extreme_Values = all_extreme_data,
      Summary = extreme_summary,
      Distribution = extreme_distribution
    ), extreme_file)
    print(paste("\n极端值详细信息已保存到:", extreme_file))
    
    # 保存最终清理后的数据
    cleaned_file <- file.path(result_path, paste0(time_window, "_kurtosis_cleaned_data.xlsx"))
    write_xlsx(current_data, cleaned_file)
    print(paste("清理后的数据已保存到:", cleaned_file))
    
    # 准备详细的峰度描述
    final_kurtosis <- kurtosis(current_data$Residual)
    kurtosis_status <- if(abs(final_kurtosis) < 2) {
      "已达到可接受范围"
    } else if(nrow(current_data) <= min_sample_size) {
      "因样本量限制停止清理"
    } else if(iteration > max_iterations) {
      "达到最大迭代次数"
    } else {
      "仍未达到理想范围"
    }
    
    kurtosis_detail <- sprintf(
      "经过%d次迭代清理，共发现%d个极端值。最终清理后峰度为%.2f（原始峰度%.2f）。状态：%s\n每次迭代情况：%s",
      iteration - 1,
      nrow(all_extreme_data),
      final_kurtosis,
      kurt_val,
      kurtosis_status,
      paste(sapply(1:max_iterations, function(i) {
        if (extreme_summary$Extreme_Values_Found[i] > 0) {
          sprintf("\n- 第%d次 (IQR倍数%.1f): 发现%d个极端值，清理后峰度%.2f，剩余样本量%d (%.1f%%)",
                 i,
                 extreme_summary$IQR_Multiplier[i],
                 extreme_summary$Extreme_Values_Found[i],
                 extreme_summary$Kurtosis_After_Removal[i],
                 extreme_summary$Remaining_Sample_Size[i],
                 100 * extreme_summary$Remaining_Sample_Size[i] / nrow(original_data))
        } else {
          ""
        }
      }), collapse = "")
    )
  } else {
    kurtosis_detail <- "峰度正常，无需处理极端值"
  }
  
  # Detailed kurtosis interpretation
  if (abs(kurt_val) < 2) {
    print("Kurtosis: Acceptable")
    kurtosis_interpretation <- "Normal distribution"
  } else if (kurt_val > 2) {
    print("Kurtosis: May have kurtosis issues (Leptokurtic)")
    kurtosis_interpretation <- "Heavy-tailed distribution (more outliers than normal)"
  } else {
    print("Kurtosis: May have kurtosis issues (Platykurtic)")
    kurtosis_interpretation <- "Light-tailed distribution (fewer outliers than normal)"
  }
  
  # Return result data frame
  result <- data.frame(
    Time_Window = time_window,
    Mean = mean(residuals),
    SD = sd(residuals),
    Skewness = skew_val,
    Kurtosis = kurt_val,
    Skewness_Acceptable = abs(skew_val) < 1,
    Kurtosis_Acceptable = abs(kurt_val) < 2,
    Mean_Acceptable = abs(mean(residuals)) < 0.1,
    Skewness_Criteria = "|value| < 1 is acceptable",
    Kurtosis_Criteria = "|value| < 2 is acceptable",
    Mean_Criteria = "|value| < 0.1 is acceptable",
    Kurtosis_Interpretation = kurtosis_interpretation,
    Kurtosis_Detail = kurtosis_detail,
    Original_Kurtosis = kurt_val,
    Final_Kurtosis = kurtosis(current_data$Residual),
    Total_Extreme_Values = nrow(all_extreme_data),
    Iterations_Used = iteration - 1,
    Final_Sample_Size = nrow(current_data),
    Removed_Data_Percentage = (nrow(all_extreme_data) / nrow(original_data)) * 100,
    stringsAsFactors = FALSE
  )
  
  return(result)
}

# Fit model and check assumptions for each time window
for (time_window in erp_columns) {
  print(paste("\nAnalyzing time window:", time_window))
  
  # Build formula
  formula <- as.formula(paste(time_window, "~ Semantic_Relatedness * Priming_Word_Type * Group + (1|Subject) + (1|Channel)"))
  
  # Fit model
  model <- lmer(formula, data = erp_data)
  
  # Check model assumptions and collect results
  current_results <- check_model_assumptions(model, time_window)
  
  results_df <- rbind(results_df, current_results)
}

# Save results directly
output_file <- file.path(result_path, "model_assumptions_results.xlsx")
write_xlsx(results_df, output_file)
print(paste("\nResults saved to:", output_file))

# Create a simple explanation file
explanations <- data.frame(
  Variable = c(
    "Time_Window", "Mean", "SD", "Skewness", "Kurtosis",
    "Skewness_Acceptable", "Kurtosis_Acceptable", "Mean_Acceptable",
    "Skewness_Criteria", "Kurtosis_Criteria", "Mean_Criteria",
    "Kurtosis_Interpretation", "Kurtosis_Detail",
    "Original_Kurtosis", "Final_Kurtosis", "Total_Extreme_Values",
    "Iterations_Used", "Final_Sample_Size", "Removed_Data_Percentage"
  ),
  Explanation = c(
    "时间窗口",
    "均值（残差的中心趋势）",
    "标准差",
    "偏度（分布的对称性）",
    "峰度（分布的尖峭程度）",
    "偏度是否可接受",
    "峰度是否可接受",
    "均值是否可接受",
    "偏度标准：|值| < 1 表示分布对称性良好",
    "峰度标准：|值| < 2 表示分布形态正常",
    "均值标准：|值| < 0.1 表示无系统性偏差",
    "峰度解释",
    "峰度清理过程详细描述",
    "原始峰度值",
    "最终清理后峰度值",
    "清理的极端值总数",
    "使用的迭代次数",
    "最终样本量",
    "被移除数据的百分比"
  )
)

# Save explanations to a separate file
explanation_file <- file.path(result_path, "variable_explanations.xlsx")
write_xlsx(explanations, explanation_file)
print(paste("\nExplanations saved to:", explanation_file))

# Print summary results
print("\n============= Model Assumption Test Summary =============")
print(results_df)
print("\n============= Variable Explanations =============")
print(explanations)







